perm filename SSS.F4[P11,LCS] blob
sn#120530 filedate 1978-03-15 generic text, type T, neo UTF8
SUBROUTINE SSS(VV,N1,A1)
DIMENSION V(50,4),A1(512),C(30,4),YP(30),J(30),NX(3),KA(14),K(9)
DIMENSION VV(50,4)
EQUIVALENCE(K1,K(1)),(K2,K(2)),(K3,K(3)),(K4,K(4)),(K5,K(5)),
1 (K6,K(6)),(K7,K(7)),(K8,K(8)),(K9,K(9))
DATA KA/1,2,2,1,1,2,1,1,0,2,1,-1,0,1/,DX/.00001/
IF(VV(1,2).EQ.0) VV(1,2)=1
DO 5 I=1,30
DO 5 L=1,2
5 V(I,L)=VV(I,L)
NX(1)=N1
698 NX(2)=NX(1)-1
DO 10 I=1,NX(1)
10 V(I,2)=(V(I,2)-1)/99.
DO 20 I=2,NX(2)
JX=I+1
JZ=I-1
YP(I)=(V(JX,1)-V(JZ,1))/(V(JX,2)-V(JZ,2))
20 IF((V(JX,1)-V(I,1))*(V(I,1)-V(JZ,1)).LE.0) YP(I)=0
DO 22 I=1,9
22 K(I)=KA(I)
KOUNT=0
21 KOUNT=KOUNT+1
V1=V(K2,1)-V(K1,1)
V2=V(K2,2)-V(K1,2)
802 IF((YP(K2)-V1/V2)*(V(K3,1)-V(K4,1)).GT.0) GO TO 30
24 Z=V(K2,K5)+(V(K1,K6)-V(K2,K6))*YP(K2)**K7
IF(YP(K2)**2.LT.DX.AND.V1**2.LT.DX) GO TO 36
IF(YP(K2)**2.LT.DX) GO TO 38
D1=V(K2,K5)-Z
806 D2=Z-V(K1,K5)
ZZ=(V(K1,K6)*D2+V(K2,K6)*D1)/(D1+D2)
808 YP(K1)=(ZZ*K9+V(K2,1)*K8-V(K1,1))/
1 (ZZ*K8+V(K2,2)*K9-V(K1,2))
GO TO 40
30 DO 32 I=5,9
32 K(I)=KA(I+5)
GO TO 24
36 YP(K1)=0
GO TO 40
38 YP(K1)=-100
IF(KOUNT.EQ.2) GO TO 39
IF(V(K2,1).GT.V(K1,1)) YP(K1)=100
GO TO 40
39 IF(V(K2,1).LT.V(K1,1)) YP(K1)=100
40 IF(KOUNT.EQ.2) GO TO 50
DO 42 I=1,2
K(I)=NX(I)
42 K(I+2)=K(I)
DO 44 I=5,9
44 K(I)=KA(I)
GO TO 21
50 NX(3)=NX(2)-1
N=1
52 N=N+1
IF(N.GT.NX(3)) GO TO 92
JX=N+1
V1=V(JX,1)-V(N,1)
V2=V(JX,2)-V(N,2)
Y1=YP(N)-YP(JX)
IF(Y1**2.LT.DX.AND.V1**2.GT.DX) GO TO 720
710 X=(V1-YP(JX)*V(JX,2)+YP(N)*V(N,2))/Y1
715 IF(X.GE.V(N,2).AND.X.LE.V(JX,2)) GO TO 52
IF(Y1**2.LT.DX.AND.V1**2.LT.DX) GO TO 52
720 DO 120 I=NX(1),JX,-1
JZ=I+1
V(JZ,2)=V(I,2)
V(JZ,1)=V(I,1)
120 YP(JZ)=YP(I)
YP(JX)=.5*V1/V2
IF(V1*(YP(N)-V1/V2).LE.0) YP(N+1)=4*YP(JX)
V(JX,2)=.5*(V(N+2,2)+V(N,2))
V(JX,1)=.5*(V(N+2,1)+V(N,1))
N=JX
DO 88 L=1,3
88 NX(L)=NX(L)+1
GO TO 52
92 DO 140 I=1,NX(2)
JX=I+1
W0=YP(I)
W1=YP(JX)
W2=V(JX,2)-V(I,2)
W3=V(JX,1)-V(I,1)
C(I,1)=(W2*(W0+W1)-2*W3)/(W0-W1)
C(I,2)=W2-C(I,1)
C(I,4)=W0*C(I,2)
140 C(I,3)=-C(I,4)+W3
730 DO 150 I=1,NX(1)
150 J(I)=511*V(I,2)+1
740 DO 160 I=1,NX(2)
L1=J(I)+1
IF(I.EQ.1) L1=1
ZZ=C(I,2)
XX=C(I,1)
L2=J(I+1)
750 DO 160 L=L1,L2
X=(FLOAT(L)-1.)/511.
IF(XX**2.LT.DX) GO TO 155
ZX=.5*SQRT(ZZ**2-4*XX*(V(I,2)-X))/XX
T1=-.5*ZZ/XX+ZX
T2=T1-2*ZX
IF(T2.GT.-DX.AND.T2.LT.(1+DX)) T1=T2
155 IF(XX**2.LT.DX) T1=-(V(I,2)-X)/ZZ
160 A1(L)=C(I,3)*T1**2+C(I,4)*T1+V(I,1)
770 END